/**
 * Implement the vector computation. 
 * The classes Vector and Matrix are provided. 
 * @author Songgang Xu (sxu at tamu dot edu)
 * @file matrix.h
 */

#ifndef XMESH_MATRIX_H
#define XMESH_MATRIX_H


namespace math
{
	template <typename value_type, unsigned dim> class vector_base 
	{
	protected:
		//std::vector<value_type> mv_vector; 		
        value_type mv_vector[dim];
		vector_base (){}
	};

	template <> class vector_base <double, 2>  
	{
	protected:
		vector_base (): x(this->mv_vector[0]), y(this->mv_vector[1]){ x = 0; y = 0;}
		double mv_vector[2];
	public:
		double & x, & y; 
	} ;

	template <> class vector_base <double, 3>  
	{
	protected:
		vector_base (): x(this->mv_vector[0]), y(this->mv_vector[1]), z(this->mv_vector[2]), r(this->mv_vector[0]), g(this->mv_vector[1]), b(this->mv_vector[2]){ x = 0; y = 0; z = 0;}
		double mv_vector[3];
	public:
		double & x, & y, & z; // for coordinate ;
		double & r, & g, & b; 
	} ;

	template <> class vector_base <double, 4>  
	{
	protected:
		vector_base (): x(this->mv_vector[0]), y(this->mv_vector[1]), 
			z(this->mv_vector[2]), w(this->mv_vector[3]){ x = 0; y = 0; z = 0; w = 0;}
		double mv_vector[4];
	public:
		double & x, & y, & z, & w; 
	} ;

	template <unsigned dim > class vector_base <double, dim>  
	{
	protected:
		//std::vector<value_type> mv_vector; 		
        double mv_vector[dim];
		vector_base (){ for (unsigned i = 0; i < dim; ++i) mv_vector[i] = 0 ;}
	};

	template <typename value_type> class vector_base <value_type, 2>  
	{
	protected:
		vector_base (): x(this->mv_vector[0]), y(this->mv_vector[1]){}
		value_type mv_vector[2];
	public:
		value_type & x, & y; 
	} ;

	template <typename value_type> class vector_base <value_type, 3> 
	{
	protected:
		vector_base () : x(this->mv_vector[0]), y(this->mv_vector[1]), z (this->mv_vector[2]), 
			r(this->mv_vector[0]), g(this->mv_vector[1]), b (this->mv_vector[2]){}
		vector_base (vector_base<value_type, 3> & _v) : x(this->mv_vector[0]), y(this->mv_vector[1]), z (this->mv_vector[2]), r(this->mv_vector[0]), g(this->mv_vector[1]), b (this->mv_vector[2])
		{
		}
		//std::vector<value_type> mv_vector; 		
        value_type mv_vector[3];
	public:
		value_type & x, & y, & z; 
		value_type & r, & g, & b; 
	} ;
	template <typename value_type> class vector_base <value_type, 4> 
	{
	protected:
		vector_base (): x(this->mv_vector[0]), y(this->mv_vector[1]), z (this->mv_vector[2]), w (this->mv_vector[3]), r(this->mv_vector[0]), g(this->mv_vector[1]), b (this->mv_vector[2]), a (this->mv_vector[3]){}
		//std::vector<value_type> mv_vector; 		//
		value_type mv_vector[4];
	public:
		value_type & x, & y, & z, & w; 
		value_type & r, & g, & b, & a; 
	} ;

	template <typename value_, unsigned dim>
		class vector : public vector_base <value_, dim>
	{
	public:
		typedef value_ value_type ;
		typedef value_type * iterator ;
		typedef const value_type * const_iterator ;
		typedef value_type &       reference ;
		typedef const value_type & const_reference ;
		typedef unsigned           size_type ;

		typedef vector<value_type, dim> vector_type ;
		typedef vector_type             point_type ;
	public:
		//itrator support 
		iterator begin() {STATIC_ASSERT (dim > 0); return this->mv_vector;}
		const_iterator begin() const {STATIC_ASSERT (dim >0) ; return this->mv_vector; }
		iterator end () {STATIC_ASSERT(dim > 0) ; return this->mv_vector + dim ;}
		const_iterator end () const {STATIC_ASSERT(dim > 0) ; return this->mv_vector + dim ;}
		typedef std::reverse_iterator<iterator>       reverse_iterator;      
		typedef std::reverse_iterator<const_iterator> const_reverse_iterator;

        reverse_iterator rbegin() { return reverse_iterator(end()); }
        const_reverse_iterator rbegin() const {
            return const_reverse_iterator(end());
        }
        const_reverse_iterator crbegin() const {
            return const_reverse_iterator(end());
        }

        reverse_iterator rend() { return reverse_iterator(begin()); }
        const_reverse_iterator rend() const {
            return const_reverse_iterator(begin());
        }
        const_reverse_iterator crend() const {
            return const_reverse_iterator(begin());
        }

		size_type size() const {return dim ;}
		bool empty () const{return dim == 0;}
		
		reference operator [] (unsigned i ) 
		{ 
			//TEST (i < dim, "Exceed index"); 
			return this->mv_vector[i]; 
		}
		const_reference operator [] (unsigned i ) const 
		{   //TEST (i < dim, "Exceed index"); 
			return this->mv_vector[i];
		}
		
		
	public:
		/**
		 * Constructor. 
		 */
		vector () : vector_base <value_type, dim> (), m_idxcount(0) {}
		vector (const value_type & _v1, const value_type & _v2) : vector_base <value_type, 2> (), m_idxcount(0)
		{
			this->mv_vector[0] = _v1; this->mv_vector[1] = _v2; 
		}
		vector (const value_type & _v1, const value_type & _v2, const value_type & _v3) : vector_base <value_type, 3> (), m_idxcount(0)  
		{
			this->mv_vector[0] = _v1; this->mv_vector[1] = _v2; this->mv_vector[2] = _v3; 
		}
		vector (const value_type & _v1, const value_type & _v2, const value_type & _v3, const value_type & _v4) : vector_base <value_type, 4> () , m_idxcount(0)
		{
			this->mv_vector[0] = _v1; this->mv_vector[1] = _v2; 
			this->mv_vector[2] = _v3; this->mv_vector[3] = _v4; 
		}

		/**
		 * Copy Constructor. 
		 * @param _v Vector instance. 
		 */
		vector (const vector <value_type, dim> & _v) : vector_base <value_type , dim>() 
		{
			for (unsigned i = 0; i < dim; ++i) this->mv_vector[i] = _v[i];
		}

		/*
		  construct from a vector . 
		template <typename element>
		vector (const std::vector <element> & _v) : vector_base <value_type , dim>() 
		{
			STATIC_ASSERT (boost::is_scalar<element>);
			for (unsigned i = 0; i < dim; ++i){this->mv_vector[i] = _v.at(i);}
		}
		*/

		/**
		 * Overload operator = 
		 */
		vector <value_type, dim > & operator = ( const vector <value_type, dim > & rhs) 
			{
				for (unsigned i = 0; i < dim; ++i) this->mv_vector[i] = rhs[i];
				return (*this); 
			}

		/**
		 * Overload operator = (value_type) for vector initialization. 
		 * @param _v Value. 
		 * @return Vector reference. 
		 */
		vector & operator = (const value_type & _v) 
			{
				TEST ( this->m_idxcount < dim, "Exceed index");
				this->m_idxcount = 0; 
				this->mv_vector[m_idxcount] = _v; 
				return (*this); 
			}

		/**
		 * Overload operator comma  for vector initialization.
		 * @param _v Value.
		 * @return Vector reference. 
		 */
		vector & operator , (const value_type & _v)
		{
			++ (this->m_idxcount ) ;
			TEST ( this->m_idxcount < dim, "Exceed index");
			this->mv_vector[m_idxcount] = _v;
			return (*this);
		}

		/**
		 * Calculate the module square of Vector. 
		 */
		value_type module2 () const {return (*this) * (*this) ;}

		/**
		 * Calculate the module of vector. 
		 */
		value_type module () const {return sqrt (module2());}

		/**
		 * Normalize the vector. 
		 */
		vector & normalize ()
		{
			value_type m = module (); 
			for (unsigned i = 0; i < dim; ++i)   this->mv_vector[i] /= m ;
			return (*this);
		}

		vector & zero ()
		{
			for (unsigned i = 0; i < dim;  ++i) this->mv_vector[i] = value_type(0);
			return (*this) ;
		}

		std::string to_string () const 
		{
			std::ostringstream out;
			out <<'[';
			//out << std::setprecision(10);
			for (unsigned i = 0; i < dim - 1; ++i) out << (*this)[i] <<',';
			if (dim > 0) out <<(*this)[dim -1]<<']';
			return out.str();
		}
	private:
		unsigned m_idxcount ;
	};

    /**
	 * Overload of ostream operator
	 * @param out Output stream.  
	 * @param _v Vector. 
	 */
	template<typename value_type, unsigned dim>
		inline std::ostream& operator<<(std::ostream & out, const vector<value_type, dim> & _v) 
	{
		out<<_v.to_string();
		return out;
	}

    /**
	 * Overload of istream operator
	 * @param in Input stream.  
	 * @param _v Vector. 
	 */
	template<typename value_type, unsigned dim>
		inline std::istream& operator >> (std::istream& in, vector<value_type, dim> & _v) 
	{
		for (unsigned i = 0; i < _v.size(); ++i) in >> _v[i];
		return in ;
	}
	
    /**
	 * Overload of operator +
	 * @param lhs Left hand operand. 
	 * @param rhs Right hand operand. 
	 */
	template<typename value_type, unsigned dim>
		inline vector <value_type, dim > operator + (const vector <value_type, dim > & lhs, const vector <value_type, dim > & rhs) 
	{
		vector <value_type, dim > result ; 
		for (unsigned i = 0; i < lhs.size(); ++i) result[i] = lhs[i] + rhs[i]; 
		return result; 
	}

	template<typename value_type, unsigned dim>
		inline vector <value_type, dim > & operator += (vector <value_type, dim > & lhs, const vector <value_type, dim > & rhs) 
	{
		for (unsigned i = 0; i < lhs.size(); ++i) lhs[i] += rhs[i];
		return lhs; 
	}

	template<typename value_type, unsigned dim>
		inline vector <value_type, dim > & operator += (vector <value_type, dim > & lhs, const value_type & rhs) 
	{
		for (unsigned i = 0; i < lhs.size(); ++i) lhs[i] += rhs;
		return lhs; 
	}

	template<typename value_type, typename vector_type, unsigned dim>
		inline vector <vector_type, dim > & operator += (vector <vector_type, dim > & lhs, const value_type & rhs) 
	{
		for (unsigned i = 0; i < lhs.size(); ++i) lhs[i] += rhs;
		return lhs; 
	}

    /**
	 * Overload of operator -
	 * @param lhs Left hand operand. 
	 * @param rhs Right hand operand. 
	 */
	template<typename value_type, unsigned dim>
		inline vector<value_type, dim> operator - (const vector<value_type, dim> & lhs, const vector<value_type, dim> & rhs) 
	{
		vector <value_type, dim > result ; 
		for (unsigned i = 0 ;i < lhs.size(); ++i) result [i] = lhs[i] - rhs[i];
		return result; 
	}

	/**
	 * Negative
	 * @param lhs. Vector. 
	 */
	template<typename value_type, unsigned dim>
		inline vector<value_type, dim>  operator - (const vector<value_type, dim> & lhs) 
	{
		vector <value_type, dim > result ; 
		for (unsigned i = 0 ;i < lhs.size(); ++i) result[i] = -lhs[i];
		return result; 
	}

	template<typename value_type, unsigned dim>
		inline vector <value_type, dim > & operator -= (vector <value_type, dim > & lhs, const vector <value_type, dim > & rhs) 
	{
		for (unsigned i = 0; i < lhs.size(); ++i) lhs[i] -= rhs[i];
		return lhs; 
	}

	template<typename value_type, unsigned dim>
		inline vector <value_type, dim > & operator -= (vector <value_type, dim > & lhs, const value_type & rhs) 
	{
		for (unsigned i = 0; i < lhs.size(); ++i) lhs[i] -= rhs;
		return lhs; 
	}

	template<typename value_type, typename vector_type, unsigned dim>
		inline vector <vector_type, dim > & operator -= (vector <vector_type, dim > & lhs, const value_type & rhs) 
	{
		for (unsigned i = 0; i < lhs.size(); ++i) lhs[i] -= rhs;
		return lhs; 
	}


    /**
	 * Overload of operator *. Dot product. 
	 * @param lhs Left hand operand. 
	 * @param rhs Right hand operand. 
	 */
	template<typename value_type, unsigned dim>
		inline value_type operator * (const vector <value_type, dim> & lhs, const vector <value_type, dim> & rhs)  
	{
		value_type result = value_type(0);
		for (unsigned i = 0 ;i < lhs.size(); ++i) result += lhs[i] * rhs[i];
		return result; 
	}

    /**
	 * Overload of operator *. vector x value. 
	 * @param lhs Left hand operand. 
	 * @param rhs Right hand operand. 
	 */
	template<typename value_type, unsigned dim>
		inline vector<value_type, dim> operator * (const vector <value_type, dim> & lhs, const value_type & rhs) 
	{
		vector<value_type, dim> result ;
		for (unsigned i = 0 ;i < lhs.size(); ++i) result[i] = lhs[i] * rhs;
		return result; 
	}

	template<typename value_type, unsigned dim>
		inline vector<value_type, dim> & operator *= (vector <value_type, dim> & lhs, const value_type & rhs)  
	{
		for (unsigned i = 0 ;i < lhs.size(); ++i) lhs[i] *= rhs;
		return lhs; 
	}

	template<typename value_type, typename vector_type, unsigned dim>
		inline vector <vector_type, dim > & operator *= (vector <vector_type, dim > & lhs, const value_type & rhs) 
	{
		for (unsigned i = 0; i < lhs.size(); ++i) lhs[i] *= rhs;
		return lhs; 
	}

    /**
	 * Overload of operator *. value x vector. 
	 * @param lhs Left hand operand. 
	 * @param rhs Right hand operand. 
	 */
	template<typename value_type, unsigned dim>
		inline vector<value_type, dim>  operator * (const value_type & lhs, const vector <value_type, dim> & rhs)
	{
		return (rhs * lhs); 
	}

    /**
	 * Overload of operator /. vector / value. 
	 * @param lhs Left hand operand. 
	 * @param rhs Right hand operand. 
	 */
	template<typename value_type, unsigned dim>
		inline vector<value_type, dim> operator / (const vector <value_type, dim> & lhs, const value_type & rhs)  
	{
		vector<value_type, dim> result ;
		for (unsigned i = 0 ;i < lhs.size(); ++i) result [i] =  lhs[i] / rhs;
		return result; 
	}

	template<typename value_type, typename vector_type, unsigned dim>
		inline vector <vector_type, dim > & operator /= (vector <vector_type, dim > & lhs, const value_type & rhs) 
	{
		for (unsigned i = 0; i < lhs.size(); ++i) lhs[i] /= rhs;
		return lhs; 
	}

    /**
	 * Overload of operator ^. Cross product of two vectors. 
	 * @param lhs Left hand operand. 3D vector.
	 * @param rhs right hand operand. 3D vector. 
	 */
	template<typename value_type>
		inline vector<value_type, 3> operator ^ (const vector<value_type, 3> & lhs, const vector<value_type, 3> & rhs) 
	{
		return vector <value_type, 3>  (lhs[1] * rhs[2] - lhs[2] * rhs[1], 
										lhs[2] * rhs[0] - lhs[0] * rhs[2],
										lhs[0] * rhs[1] - lhs[1] * rhs[0] );
	}

    /**
	 * Overload of operator ^. Cross product of two vectors. 
	 * @param lhs Left hand operand. 2D vector.
	 * @param rhs right hand operand. 2D vector. 
	 */
	template<typename value_type>
		inline value_type operator ^ (const vector<value_type, 2> & lhs, const vector<value_type, 2> & rhs) 
	{
		return (lhs[0] * rhs[1] - lhs[1] * rhs[0] );
	}

	/**
	 * Overload comparison operator <.
	 * We order two vectors having the same dimension in the alphabetic order. 
	 */
	template<typename value_type, unsigned dim >
		inline bool operator < (const vector<value_type, dim> & lhs, const vector<value_type, dim> & rhs)
	{
		for (unsigned i = 0; i < dim; ++i){
			if (lhs[i] < rhs[i]) return true; 
			if (lhs[i] > rhs[i]) return false; 
		}
		return false; 
	}

    /**
	 * Overload comparison operator >.
	 * We order two vectors having the same dimension in the alphabetic order. 
	 */
	template<typename value_type, unsigned dim >
		inline bool operator > (const vector<value_type, dim> & lhs, const vector<value_type, dim> & rhs)
	{return rhs < lhs; }

    /**
	 * Overload comparison operator <=.
	 * We order two vectors having the same dimension in the alphabetic order. 
	 */
	template<typename value_type, unsigned dim >
		inline bool operator <= (const vector<value_type, dim> & lhs, const vector<value_type, dim> & rhs)
	{return !(rhs < lhs); }


    /**
	 * Overload comparison operator >=.
	 * We order two vectors having the same dimension in the alphabetic order. 
	 */
	template<typename value_type, unsigned dim >
		inline bool operator >= (const vector<value_type, dim> & lhs, const vector<value_type, dim> & rhs)
	{return !(lhs < rhs); }

	/**
	 * Overload comparison operator ==.
	 * We order two vectors having the same dimension in the alphabetic order. 
	 */
	template<typename value_type, unsigned dim >
		inline bool operator == (const vector<value_type, dim> & lhs, const vector<value_type, dim> & rhs)
	{
		for (unsigned i = 0; i < dim; ++i) if (lhs[i] != rhs[i]) return false; 
		return true; 
	}

	/**
	 * Overload comparison operator !=.
	 * We order two vectors having the same dimension in the alphabetic order. 
	 */
	template<typename value_type, unsigned dim >
		inline bool operator != (const vector<value_type, dim> & lhs, const vector<value_type, dim> & rhs)
	{
		return !(lhs == rhs);
	}

	/**
	 * functions. 
	 */
	template<typename value_type, unsigned dim>
		inline vector <value_type, dim > abs (const vector <value_type, dim > & lhs)
	{
		vector <value_type, dim> result ;
		for (unsigned i = 0; i < lhs.size(); ++i) result[i] = std::abs(lhs[i]);
		return result; 
	}


	template<typename value_type, unsigned dim>
		inline vector <value_type, dim > min (const vector <value_type, dim > & lhs, const vector <value_type, dim > & rhs)
	{
		vector <value_type, dim> result ;
		for (unsigned i = 0; i < lhs.size(); ++i) result[i] = std::min(lhs[i], rhs[i]);
		return result; 
	}

	template<typename value_type, unsigned dim>
		inline vector <value_type, dim > max (const vector <value_type, dim > & lhs, const vector <value_type, dim > & rhs)
	{
		vector <value_type, dim> result ;
		for (unsigned i = 0; i < lhs.size(); ++i) result[i] = std::max(lhs[i], rhs[i]);
		return result; 
	}
    
	template<typename value_type, unsigned dim>
		inline value_type distance (const vector<value_type, dim> & lhs, const vector<value_type, dim> & rhs)
	{
		return (lhs-rhs).module();
	}

	template<typename value_type, unsigned dim>
		inline value_type distance2 (const vector<value_type, dim> & lhs, const vector<value_type, dim> & rhs)
	{
		return (lhs-rhs).module2();
	}

	// Compute lhs^ rhs
	template<typename value_type>
		inline value_type orient(const vector<value_type, 2> &lhs, const vector<value_type, 2> &rhs)
	{
		return (lhs^rhs);
    }

	// Compute a . (b x c)
	template<typename value_type>
		inline value_type orient(const vector<value_type, 3> & a, const vector<value_type, 3> & b, const vector<value_type, 3> & c) 
	{
		return
			(a.x * b.y * c.z + a.y * b.z * c.x + a.z * b.x * c.y) -
			(a.x * b.z * c.y + a.y * b.x * c.z + a.z * b.y * c.x);
    }

	// select the index has the abs maximum value. 
	template<typename value_type, unsigned dim>
		unsigned max_index (const vector <value_type, dim> & _v)
	{
		STATIC_ASSERT(dim > 0);
		vector<value_type, dim> v = abs(_v); 
		unsigned coord = 0; 
		for (unsigned i = 1; i < dim; ++i)
		{
			if (v[i] > v[coord]) coord = i;
		}
		return coord; 
	}

	template<typename value_type, unsigned row, unsigned col> class matrix  ;
	/**
	 * Matrix base
	 */
	template<typename value_type, unsigned row, unsigned col> 
		class matrix_base 
	{
	protected:
		typedef vector <value_type, col> row_vector_type ;
		vector< row_vector_type, row> mv_matrix;
	};

	template<typename value_type>
		class matrix_base <value_type, 2, 2>
	{
	public:
		value_type det () const { return  mv_matrix[0][0] * mv_matrix[1][1] - mv_matrix[0][1] * mv_matrix[1][0]; }
		matrix <value_type, 2, 2> inv ()
		{
			matrix <value_type , 2, 2> invm ;
			const vector< row_vector_type, 2> & m = mv_matrix ;
			value_type rdet = value_type(1.0) / this->det();
			invm[0][0] = rdet * m[1][1]; invm[0][1] = - rdet * m[1][0] ;
			invm[1][0] = -rdet * m[0][1]; invm[1][1] = rdet * m[0][0] ;
			return invm; 
		}
	protected:
		typedef vector <value_type, 2> row_vector_type ;
		vector< row_vector_type, 2> mv_matrix;
	};

	template<typename value_type>
		class matrix_base <value_type, 3, 3>
	{
	public:
		value_type det () const { 
			const vector< row_vector_type, 3> & m = mv_matrix ;
			return ( m[0][0] *(m[2][2]*m[1][1] - m[2][1]* m[1][2]) 
					 - m[1][0] * (m[2][2]*m[0][1] - m[2][1]* m[0][2])
					 + m[2][0] * (m[1][2]*m[0][1] - m[1][1]* m[0][2]));
		}
		matrix <value_type, 3, 3> inv () const 
		{
			matrix <value_type , 3, 3> invm ;
			const vector< row_vector_type, 3> & m = mv_matrix ;
			value_type rdet = value_type(1.0) / this->det();
			invm[0][0] =   rdet * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) ;
			invm[0][1] = - rdet * (m[0][1] * m[2][2] - m[0][2] * m[2][1]) ;
			invm[0][2] =   rdet * (m[0][1] * m[1][2] - m[0][2] * m[1][1]) ;

			invm[1][0] = - rdet * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) ;
			invm[1][1] =   rdet * (m[0][0] * m[2][2] - m[0][2] * m[2][0]) ;
			invm[1][2] = - rdet * (m[0][0] * m[1][2] - m[0][2] * m[1][0]) ;

			invm[2][0] =   rdet * (m[1][0] * m[2][1] - m[1][1] * m[2][0]) ;
			invm[2][1] = - rdet * (m[0][0] * m[2][1] - m[0][1] * m[2][0]) ;
			invm[2][2] =   rdet * (m[0][0] * m[1][1] - m[0][1] * m[1][0]) ;
			return invm; 
		}
	protected:
		typedef vector <value_type, 3> row_vector_type ;
		vector< row_vector_type, 3> mv_matrix;
	};

	/**
	 * Matrix class 
	 */
	template<typename value_type, unsigned row, unsigned col>
		class matrix : public matrix_base <value_type, row, col>
	{
	protected:
		typedef vector <value_type, col> row_vector_type ; 
		typedef vector <value_type, row> col_vector_type ; 
		typedef matrix <value_type, row, col> matrix_type ;
		typedef matrix <value_type, col, row> trans_matrix_type ;
	public:
		matrix () : m_idxcount(0){}
		matrix (const matrix_type & _m) : m_idxcount(0){ this->mv_matrix = _m.mv_matrix ;}
		matrix_type & operator = (const matrix_type & _m) {this->mv_matrix = _m.mv_matrix; return (*this);}
		matrix_type & operator = (const value_type & _v) 
			{
				STATIC_ASSERT (row != 0 || col !=0);
				this->m_idxcount = 0;
				this->mv_matrix[0][0] = _v; 
				return (*this); 
			}
		matrix_type & operator , (const value_type & _v)
		{
			++ this->m_idxcount ;
			TEST (this->m_idxcount < row * col, "Exceed Index") ;			
			unsigned crow = m_idxcount / col ;
			unsigned ccol = m_idxcount % col ;
			this->mv_matrix[crow][ccol] = _v; 
			return (*this);
		}

		unsigned row_size() const {return row; }
		unsigned col_size() const {return col; }

		row_vector_type & operator [] (unsigned _i ) {return this->mv_matrix[_i];}
		const row_vector_type & operator [] (unsigned _i ) const {return this->mv_matrix[_i];}
		value_type & operator () (unsigned _r, unsigned _c) {return this->mv_matrix[_r][_c];} 
		const value_type & operator () (unsigned _r, unsigned _c) const {return this->mv_matrix[_r][_c];}

		matrix_type & zero () {for (unsigned i = 0; i < row; ++i) this->mv_matrix[i].zero(); return (*this);}
		bool is_square () const {return row == col ;}

		trans_matrix_type  transpose () const 
		{
			trans_matrix_type  trans; 
			for (unsigned i = 0; i < row; ++i)
				for (unsigned j = 0; j < col; ++j)
					trans (j, i) = this->mv_matrix[i][j] ;
			return trans; 
		}
		
		friend std::ostream & operator << (std::ostream & ostr, const matrix_type & _m)
		{
			for (unsigned i = 0; i < row; ++i) ostr<< _m[i]<<'\n';
			return ostr; 
		}

		friend std::ostream & operator >> (std::ostream & is, const matrix_type & _m)
		{
			for (unsigned i = 0; i < _m.m_row; ++i) is >> _m[i]; 
			return is;
		}

		friend matrix_type operator + (const matrix_type & lhs, const matrix_type & rhs)
		{
			matrix_type result; 
			result.mv_matrix = lhs.mv_matrix + rhs.mv_matrix; 
			return result ;
		}

		friend matrix_type & operator += (matrix_type & lhs, const matrix_type & rhs)
		{
			lhs.mv_matrix += rhs.mv_matrix; 
			return lhs ;
		}
		
		friend matrix_type operator - (const matrix_type & lhs, const matrix_type & rhs)
		{
			matrix_type result; 
			result.mv_matrix = lhs.mv_matrix - rhs.mv_matrix; 
			return result ;
		}

		friend matrix_type & operator -= (matrix_type & lhs, const matrix_type & rhs)
		{
			lhs.mv_matrix -= rhs.mv_matrix; 
			return lhs ;
		}

		friend matrix_type operator * (const matrix_type & lhs, const value_type & _v) 
		{
			matrix_type result (lhs);
			result *= _v; 
			return result ;
		}

		friend row_vector_type operator * (const matrix_type & lhs, const col_vector_type & _v)
		{
			row_vector_type result;
			for (unsigned i = 0 ; i < lhs.row_size(); ++i)
			{
				result[i] = lhs.mv_matrix[i] * _v;
			}

			return result ;

		}
		friend matrix_type operator * ( const value_type & _v, const matrix_type & rhs) 
		{
			return (rhs * _v);
		}

		friend matrix_type & operator *= (matrix_type & lhs, const value_type & _v)
		{
			lhs.mv_matrix *= _v; 
			return lhs ;
		}



		friend matrix_type operator / (const matrix_type & lhs, const value_type & _v) 
		{
			matrix_type result;
			result.mv_matrix /= _v; 
			return result ;
		}

		friend matrix_type & operator /= (matrix_type & lhs, const value_type & _v)
		{
			lhs.mv_matrix /= _v; 
			return lhs ;
		}
	private:
		unsigned m_idxcount ; 
	};
	/*
	 * Matrix multiplication. 
	 */
	template <typename value_type, unsigned row, unsigned col, unsigned col2>
		matrix <value_type, row, col2> operator * (const matrix <value_type, row, col> & lhs, const matrix <value_type, col, col2> & rhs) 
	{
		matrix <value_type, row, col2> result;
		for (unsigned i = 0; i < row; ++i)
		{
			for (unsigned j = 0; j < col2; ++j)
			{
				result(i,j) = value_type(0); 
				for (unsigned k = 0; k < col ; ++k)
				{
					result(i,j) += lhs(i,k ) * rhs(k,j);
				}
			}
		}
		return result ;
	}



    /**
     * Define Dynamic vector.  
     */
	template <typename value_>
		class dynamic_vector
	{
	protected:
		typedef value_ value_type; 
		typedef dynamic_vector<value_type> vector_type ;
	public:
		typedef dynamic_vector<value_type> point_type ;
	public:
		/**
		 * Constructor. 
		 */
		dynamic_vector ():m_size(0) {}

		/**
		 * Copy Constructor. 
		 * @param _v dynamic_vector instance. 
		 */
		dynamic_vector (const dynamic_vector & _v): m_size(_v.m_vector.size()), m_vector (_v.m_vector)  {}

		/**
		 * Constructor. 
		 * @param _s Size of vector. 
		 */
		dynamic_vector (unsigned _s): m_size(_s) {m_vector.resize(m_size); }

		/**
		 * Constructor for 2D vector. 
		 *
		 */
		dynamic_vector (const value_type & _v1, const value_type & _v2) : m_size(2)
		{
			m_vector.resize(m_size);
			m_vector[0] = _v1 ;  m_vector[1] = _v2 ; 
		}

		/**
		 * Constructor for 3D vector. 
		 *
		 */
		dynamic_vector (const value_type & _v1, const value_type & _v2, const value_type & _v3) : m_size(3)
		{
			m_vector.resize(m_size);
			m_vector[0] = _v1 ;  m_vector[1] = _v2 ; m_vector[2] = _v3 ;
		}

		/**
		 * Constructor for nD vector. 
		 */
		dynamic_vector ( unsigned _n , const value_type * _v) : m_size (_n)
		{
			m_vector.resize(m_size) ;
			for (unsigned i = 0; i < m_size; ++i)
			{
				m_vector[i] = _v[i]; 
			}
		}

		/**
		 * Get vector size. 
		 */
		unsigned size() const { return this->m_size;}

		/**
		 * Resize the vector.
		 * @param _s
		 */
		void resize (unsigned _s)
		{
			this->m_vector.resize(_s); 
			this->m_size = _s ;
		} 
	 
		/**
		 * Clear the vertor. 
		 */
		void clear ()
		{
			this->m_vector.clear();
			this->m_size = 0; 
		}

		/**
		 * Reference of vector. 
		 * @param _i Index of element. 
		 */
		value_type & operator [] (unsigned _i) {return this->m_vector[_i];}

		/**
		 * Value of element in vector. 
		 * @param _i Index of element. 
		 */
		const value_type & operator [] (unsigned _i) const {return this->m_vector[_i];}

		/**
		 * Overload operator = (ValueT) for vector initialization. 
		 * @param _v Value. 
		 * @return dynamic_vector reference. 
		 */
		dynamic_vector & operator = (const value_type & _v) 
			{
//		this->m_idxcount = 0; 
				this->m_vector.clear(); 
				this->m_vector.push_back( _v );
				this->m_size = 1 ;
				return (*this); 
			}

		/**
		 * Overload operator comma  for vector initialization.
		 * @param _v Value.
		 * @return dynamic_vector reference. 
		 */
		dynamic_vector & operator , (const value_type & _v)
		{
//		++this->m_idxcount ;
//		if (this->m_idxcount >= this->m_size) throw ("dynamic_vector initialization exceeds bound. ");
//		this->m_vector[m_idxcount] = _v; 
			this->m_vector.push_back( _v );
			++ (this->m_size ) ;
			return (*this);
		}

		/**
		 * Append a value into the vector.
		 * @param _v Value. 
		 */
		void append (const value_type & _v)
		{
			this->m_vector.push_back(_v);
			++ (this->m_size) ;
		}

		/**
		 * Make every component to be non-negative. 
		 * @return None.
		 */
		void positive ()
		{
			for (unsigned i = 0; i < m_size; ++i)
			{
				m_vector[i] = abs (m_vector[i]);
			}
		}
	public:

		/**
		 * Calculate the module square of dynamic_vector. 
		 */
		value_type square_module ()
		{
			value_type sum = value_type(0);
			for (unsigned i = 0; i < m_size ; ++i)
			{
				sum += this->m_vector[i] * this->m_vector[i]; 
			}
			return sum ;
		}

		/**
		 * Calculate the module of vector. 
		 */
		value_type module ()
		{
			return sqrt (square_module());
		}

		/**
		 * Normalize the vector. 
		 */
		dynamic_vector & normalize ()
		{
			value_type m = module (); 
			for (unsigned i = 0; i < m_size; ++i)
			{
				this->m_vector[i] /= m ;
			}
			return (*this);
		}
	


		/**
		 * Overload of operator *
		 * @param lhs Left hand operand. 
		 * @param rhs Right hand operand. 
		 */
		friend inline value_type operator * (const dynamic_vector & lhs, const dynamic_vector & rhs)  
		{
			if (!(lhs.size() == rhs.size())) throw ("The size of vectors doesn't match. ");

			value_type result = value_type(0);
			for (unsigned i = 0 ;i < lhs.size(); ++i)
			{
				result += lhs[i] * rhs[i];
			}
			return result; 
		}

		/**
		 * Overload of operator +=
		 * @param _v dynamic_vector. 
		 * @return incremented dynamic_vector. 
		 */
		inline dynamic_vector & operator += (const dynamic_vector & _v)
		{
			if (!(this->size() == _v.size())) throw ("The size of vectors doesn't match. ");
			for (unsigned i = 0 ;i < this->size(); ++i)
			{
				this->m_vector[i] += _v[i];
			}
			return (*this); 
		}

		/**
		 * Overload of operator +
		 * @param lhs Left hand operand. 
		 * @param rhs Right hand operand. 
		 */
		friend inline dynamic_vector operator + (const dynamic_vector & lhs, const dynamic_vector & rhs) 
		{
			if (!(lhs.size() == rhs.size())) throw ("The size of vectors doesn't match. ");
			dynamic_vector result (lhs);
			result += rhs; 
			return result; 
		}

		/**
		 * Overload of operator -
		 * @param lhs Left hand operand. 
		 * @param rhs Right hand operand. 
		 */
		friend inline dynamic_vector operator - (const dynamic_vector & lhs, const dynamic_vector & rhs) 
		{
			if (!(lhs.size() == rhs.size())) throw ("The size of vectors doesn't match. ");
			dynamic_vector result (lhs.size());
			for (unsigned i = 0 ;i < lhs.size(); ++i)
			{
				result [i] = lhs[i] - rhs[i];
			}
			return result; 
		}

		/**
		 * Overload of operator *
		 * @param lhs Left hand operand. 
		 * @param _v Scale value. 
		 */
		friend inline dynamic_vector operator * (const dynamic_vector & lhs, const value_type & _v ) 
		{
			dynamic_vector result (lhs.size());
			for (unsigned i = 0 ;i < lhs.size(); ++i)
			{
				result [i] = lhs[i] * _v;
			}
			return result; 
		}

		/**
		 * Overload of operator *
		 * @param _v Scale value. 
		 * @param rhs Left hand operand. 
		 */
		friend inline dynamic_vector operator * ( const value_type & _v, const dynamic_vector & rhs )
		{
			dynamic_vector result (rhs.size());
			for (unsigned i = 0 ;i < rhs.size(); ++i)
			{
				result [i] = rhs[i] * _v;
			}
			return result; 
		}

		/**
		 * Overload of operator /=
		 * @param _v dynamic_vector. 
		 * @return incremented dynamic_vector. 
		 */
		inline dynamic_vector & operator /= (const value_type & _v)
		{
			for (unsigned i = 0 ;i < this->size(); ++i)
			{
				this->m_vector[i] /= _v;
			}
			return (*this); 
		}

		/**
		 * Overload of operator /
		 * @param lhs Left hand operand. 
		 * @param _v Scale value. 
		 */
		friend inline dynamic_vector operator / (const dynamic_vector & lhs, const value_type & _v ) 
		{
			dynamic_vector result (lhs);
			result /= _v; 
			return result; 
		}

		/**
		 * Overload of operator ^. Cross product of two vectors. 
		 * @param lhs Left hand operand. 3D vector.
		 * @param rhs right hand operand. 3D vector. 
		 */
		friend inline dynamic_vector operator ^ (const dynamic_vector & lhs, const dynamic_vector & rhs) 
		{
			if (!(lhs.size() == 3 && rhs.size() == 3)) throw ("Cross product is computed from two 3D vectors. ");
			return dynamic_vector (lhs[1] * rhs[2] - lhs[2] * rhs[1], 
						   lhs[2] * rhs[0] - lhs[0] * rhs[2],
						   lhs[0] * rhs[1] - lhs[1] * rhs[0] );
		}

		/**
		 * Overload of ostream operator
		 * @param out Output stream.  
		 * @param _v dynamic_vector. 
		 */
		friend inline std::ostream& operator<<(std::ostream & out, const dynamic_vector & _v) 
		{
			out <<'[';
			for (unsigned i = 0; i < _v.m_size - 1; ++i) out << _v.m_vector[i] <<',';
			if (_v.m_size > 0) out <<_v.m_vector[_v.m_size -1]<<']';
			return out;
		}

		/**
		 * Overload of istream operator
		 * @param in Input stream.  
		 * @param _v dynamic_vector. 
		 */
		friend inline std::istream& operator>>(std::istream& in, dynamic_vector& _v) 
		{
			for (unsigned i = 0; i < _v.size(); ++i)
			{
				in >> _v.m_vector[i]; 
			}
			return in ;
		}

		/**
		 * Overload comparison operator <.
		 * We order two vectors having the same dimension in the alphabetic order. 
		 */
		friend inline bool operator < (const dynamic_vector & lhs, const dynamic_vector & rhs)
		{
			for (unsigned i = 0; i < lhs.size(); ++i)
			{
				if (lhs[i] < rhs[i]) return true; 
				if (lhs[i] > rhs[i]) return false; 
			}
			return false; 
		}
	
		/**
		 * Overload comparison operator >.
		 * We order two vectors having the same dimension in the alphabetic order. 
		 */
		friend inline bool operator > (const dynamic_vector & lhs, const dynamic_vector & rhs)
		{
			return  (rhs < lhs) ;
		}

		/**
		 * Overload of operator == 
		 * @param lhs Left hand operand. 
		 * @param rhs Right hand operand. 
		 */
		friend inline bool operator == (const dynamic_vector & lhs, const dynamic_vector & rhs)  
		{
			return (lhs.m_vector == rhs.m_vector);
		}

		/**
		 * Overload of operator != 
		 * @param lhs Left hand operand. 
		 * @param rhs Right hand operand. 
		 */
		friend inline bool operator != (const dynamic_vector & lhs, const dynamic_vector & rhs)  
		{
			return (!(lhs == rhs));
		}

		/**
		 * Overload comparison operator <=.
		 * We order two vectors having the same dimension in the alphabetic order. 
		 */
		friend inline bool operator <= (const dynamic_vector & lhs, const dynamic_vector & rhs)
		{
			return  !(rhs < lhs) ;
		}

		/**
		 * Overload comparison operator >=.
		 * We order two vectors having the same dimension in the alphabetic order. 
		 */
		friend inline bool operator >= (const dynamic_vector & lhs, const dynamic_vector & rhs)
		{
			return  !(lhs < rhs) ;
		}
	private:
		unsigned m_size; 
		std::vector<value_type> m_vector; 
	
		unsigned m_idxcount; 
	};



/**
 * Define a 2D matrix class. 
 */
	template <typename value_>
		class dynamic_matrix
	{
	protected:
		/** 
		 * Type name for value 
		 */
		typedef value_ value_type ;
		typedef dynamic_vector<value_type> vector_type; 

	public:

		/**
		 * Constructs a matrix class. 
		 */
		dynamic_matrix (): m_row (0), m_col(0){}

		/**
		 * Constructs a matrix class of r x c.
		 * @param _r row of the matrix. 
		 * @param _c col of the matrix. 
		 */
		dynamic_matrix (unsigned _r, unsigned _c) : m_row(_r), m_col (_c) 
		{
			allocate_memory();
		}

		/**
		 * Constructs a matrix class of r x c with the initial value v.
		 * @param _r row of the matrix. 
		 * @param _c col of the matrix. 
		 * @param _v initial value. 
		 */
		dynamic_matrix (int _r, int _c,  const value_type  & _v) : m_row(_r), m_col (_c) 
		{
			allocate_memory (this->m_pElements, m_row, m_col, _v);
		}

		/**
		 * Constructs a matrix from a vector. The vector is viewed to be column vector. 
		 */
		dynamic_matrix (const vector_type & _v): m_row (_v.size()), m_col (1)
		{
			allocate_memory();
			for (unsigned i = 0; i < m_row; ++i)
			{
				m_pElements[i][0] = _v[i];
			}
		}

		/**
		 * Copy Constructor.
		 * @param _m dynamic_matrix Object. 
		 */
		dynamic_matrix (const dynamic_matrix & _m)
		{
			m_row = _m.get_row();
			m_col = _m.get_col();
			allocate_memory(); 
			for (unsigned i = 0; i < m_row; ++i)
			{
				for (unsigned j = 0; j < m_col; ++j)
				{
					m_pElements[i][j] = _m.m_pElements[i][j];
				}
			}
		}

		/**
		 * Overload of operator = 
		 * @param _m dynamic_matrix Object. 
		 */

		dynamic_matrix & operator = (const dynamic_matrix & _m)
			{
				//if (this->m_pElements != NULL) deallocate_memory();
				m_row = _m.get_row();
				m_col = _m.get_col();
				allocate_memory(); 
				for (unsigned i = 0; i < m_row; ++i)
				{
					for (unsigned j = 0; j < m_col; ++j)
					{
						m_pElements[i][j] = _m.m_pElements[i][j];
					}
				}
				return (*this) ;
			}

		/**
		 * Destructor. 
		 */
		virtual ~dynamic_matrix () {deallocate_memory(); }


		/**
		 * Resize matrix. 
		 * @param _row New row value.
		 * @param _col New col value. 
		 * @return None. 
		 */
		void resize (unsigned _row, unsigned _col)
		{
			/**Fix me here*/
			this->m_row = _row;
			this->m_col = _col; 
			this->m_pElements.resize(_row);
			for (unsigned i = 0; i < _row; ++i)
			{
				m_pElements[i].resize(_col);
			}
		}

		/**
		 * Overload operator [].
		 * @param _i Index of the first dimension
		 * @return Reference of row vector. 
		 */
		dynamic_vector<value_type> & operator [] (unsigned _i)
		{
			return m_pElements[_i];
		}

		/**
		 * Overload operator [].
		 * @param _i Index of the first dimension
		 * @return Constant reference of row vector. 
		 */
		const dynamic_vector<value_type> & operator [] (unsigned _i) const 
		{
			return m_pElements[_i];
		} 


		/**
		 * Overload operator () to fetch element as L-Value. 
		 * @param _r Row value.
		 * @param _c Col value.
		 * @return Reference of element. 
		 */
		value_type & operator () (unsigned _r, unsigned _c)
		{
			return m_pElements[_r][_c];
		} 
		/**
		 * Overload operator () to fetch element. 
		 * @param _r Row value.
		 * @param _c Col value.
		 * @return Constant reference of element. 
		 */
		const value_type & operator () (unsigned _r, unsigned _c) const
		{
			return m_pElements[_r][_c];
		}
		/**
		 * Get the row of matrix. 
		 * @return the row count of matrix. 
		 */
		unsigned get_row () const {return this->m_row; }

		/**
		 * Get the col of matrix. 
		 * @return the col count of matrix. 
		 */
		unsigned get_col () const {return this->m_col; }

		/**
		 * Set the row of matrix. 
		 * @param _r Row value. 
		 */
		void set_row (unsigned _r ) {this->m_row = _r; }

		/**
		 * Set the col of matrix.
		 * @param _c Col value. 
		 */
		void set_col (unsigned _c ) {this->m_col = _c; }
	
		//unsigned get_dim1 () const {return this->m_row; }
		//unsigned get_dim2 () const {return this->m_col; }

		/**
		 * Check if the matrix is empty. 
		 * @return if the matrix is empty. 
		 */
		bool is_empty() const 
		{
			return (m_pElements == NULL && m_row == 0 && m_col == 0) ;
		}

		/**
		 * Check if the matrix is square.
		 * @return if the matrix is square. 
		 */
		bool is_square () const 
		{
			return (m_row == m_col);
		}

		/**
		 * Check if two matrices have the same dimension. 
		 * @return if two matrices have the same dimension.
		 */
		friend bool is_dim_match(const dynamic_matrix & _m1, const dynamic_matrix & _m2)
		{
			return (_m1.get_row () == _m2.get_row() && _m1.get_col () == _m2.get_col());
		}

	public:
		/**
		 * Overload plus operator. 
		 * @param lhs Left hand operand.
		 * @param rhs Right hand operand. 
		 * @return the sum of two matrix. 
		 */
		friend dynamic_matrix operator + (const dynamic_matrix & lhs, const dynamic_matrix & rhs)
		{
			if (! (lhs.get_row() == rhs.get_row() && lhs.get_col() == rhs.get_col() ))
			{
				throw ("Matrices dimensions do not match.");
			}
			unsigned _row = lhs.get_row(); 
			unsigned _col = lhs.get_col ();
			dynamic_matrix result(_row, _col); 
			for (int i = 0; i < _row; ++i)
			{
				for (int j = 0; j < _col; ++j)
				{
					result (i,j ) = lhs(i,j ) + rhs(i,j);
				}
			}
			return result ;
		}

		/**
		 * Overload minus operator. 
		 * @param lhs Left hand operand.
		 * @param rhs Right hand operand. 
		 * @return the difference of two matrix. 
		 */
		friend dynamic_matrix operator - (const dynamic_matrix & lhs, const dynamic_matrix & rhs) 
		{
			if (! (lhs.get_row() == rhs.get_row() && lhs.get_col() == rhs.get_col() ))
			{
				throw ("Matrices dimensions do not match.");
			}
			unsigned _row = lhs.get_row(); 
			unsigned _col = lhs.get_col ();
			dynamic_matrix result(_row, _col); 
			for (int i = 0; i < _row; ++i)
			{
				for (int j = 0; j < _col; ++j)
				{
					result (i,j ) = lhs(i,j ) - rhs(i,j);
				}
			}
			return result ;
		}

		/**
		 * Overload times operator. 
		 * @param lhs Left hand operand.
		 * @param rhs Scale value. 
		 * @return the product of the matrix and the scale value. 
		 */
		friend dynamic_matrix operator * (const dynamic_matrix & lhs, value_type & _v) 
		{
			unsigned _row = lhs.get_row(); 
			unsigned _col = lhs.get_col ();
			dynamic_matrix result(_row, _col); 
			for (int i = 0; i < _row; ++i)
			{
				for (int j = 0; j < _col; ++j)
				{
					result (i,j ) = lhs(i,j ) * _v;
				}
			}

			return result ;
		}

		/**
		 * Overload times operator. 
		 * @param lhs Scale value. 
		 * @param rhs Right hand operand. 
		 * @return the product of the matrix and the scale value. 
		 */
		friend dynamic_matrix operator * (value_type & _v, const dynamic_matrix & rhs) 
		{
			unsigned _row = rhs.get_row(); 
			unsigned _col = rhs.get_col ();
			dynamic_matrix result(_row, _col); 
			for (int i = 0; i < _row; ++i)
			{
				for (int j = 0; j < _col; ++j)
				{
					result (i,j ) = rhs(i,j ) * _v;
				}
			}
			return result ;
		}


		/**
		 * Overload times operator. 
		 * @param lhs Matrix.
		 * @param rhs dynamic_vector. 
		 * @return the product of the matrix and the vector. 
		 */
		friend dynamic_matrix operator * (const dynamic_matrix & lhs, const vector_type & _v)
		{
			unsigned _row = lhs.get_row(); 
			unsigned _col = lhs.get_col ();
			if ( !(lhs.get_col() == _v.size())) 
			{
				throw ("Dimensions of matrix and vector do not match. ");
			}
			dynamic_matrix vector_matrix (_v);

			return (lhs * vector_matrix);
		}

		/**
		 * Overload times operator. 
		 * @param lhs Left hand operand.
		 * @param rhs Right hand operand. 
		 * @return the product of two matrix. 
		 */
		friend dynamic_matrix operator * (const dynamic_matrix & lhs, const dynamic_matrix & rhs) 
		{
			if (! (lhs.get_col() == rhs.get_row() ))
			{
				throw ("Matrices dimensions do not match.");
			}

			unsigned _row = lhs.get_row(); 
			unsigned _col = rhs.get_col ();
			unsigned _k = lhs.get_col ();
			dynamic_matrix result(_row, _col); 
			for (unsigned i = 0; i < _row; ++i)
			{
				for (unsigned j = 0; j < _col; ++j)
				{
					result(i,j) = 0; 
					for (unsigned k = 0; k < _k ; ++k)
					{
						result(i,j) += lhs(i,k ) * rhs(k,j);
					}
				}
			}
			return result ;
		}
	
		/**
		 * Transpose the matrix.
		 * @param Matrix to be transposed.
		 * @return Transposed matrix. 
		 */
		friend dynamic_matrix transpose (const dynamic_matrix & _m) 
		{
			unsigned _col = _m.get_row(); 
			unsigned _row = _m.get_col ();
			dynamic_matrix result(_row, _col); 
			for (int i = 0; i < _row; ++i)
			{
				for (int j = 0; j < _col; ++j)
				{
					result (j, i) = _m(i, j) ;
				}
			}
			return result ;
		}


	public: 

		/**
		 * Output stream operator. 
		 * @param os Output stream.
		 * @param _m Matrix. 
		 * @return Output stream reference. 
		 */
		friend std::ostream & operator << (std::ostream & os, const dynamic_matrix & _m)
		{
			for (unsigned i = 0; i < _m.m_row; ++i)
			{
				for (unsigned j = 0; j < _m.m_col; ++j)
				{
					os << _m(i,j)<<' ';
				}
				os<<std::endl;
			}
			os<<std::endl;
			return os ;
		}


		/**
		 * Output stream operator. 
		 * @param is Input stream.
		 * @param _m Matrix. 
		 * @return Output stream reference. 
		 */
		friend std::ostream & operator >> (std::ostream & is, const dynamic_matrix & _m)
		{
			for (unsigned i = 0; i < _m.m_row; ++i)
			{
				is >> _m[i]; 
			}
			return is;
		}

		/**
		 * Overload operator = (value_type) for matrix initialization. 
		 * @param _v Value. 
		 * @return Matrix reference. 
		 */
		dynamic_matrix & operator = (const value_type & _v) 
			{
				this->m_idxcount = 0;
				m_pElements[0][0] = _v; 
				return (*this); 
			}

		/**
		 * Overload operator comma  for matrix initialization.
		 * @param _v Value.
		 * @return Matrix reference. 
		 */
		dynamic_matrix & operator , (const value_type & _v)
		{
			++ this->m_idxcount ;
			if (this->m_idxcount >= this->m_row * this->m_col) 
			{
				throw ("Initialization exceeds the matrix dimension bounds."); 
				return (*this);
			}
			unsigned row = m_idxcount / m_col ;
			unsigned col = m_idxcount % m_col ;
			m_pElements[row][col] = _v; 

			return (*this);
		}

		/**
		 * Determinant of 3x3 matrix. Dimension of other matrix is not applicable. The lu_decomposor should be used instead. 
		 */
		value_type det2 () const
		{
			const 	dynamic_vector<dynamic_vector<value_type> > & m = m_pElements ;
			value_type det = m[0][0] * m[1][1] - m[0][1] * m[1][0]; 
			return det ;
		}

		/**
		 * Determinant of 3x3 matrix. Dimension of other matrix is not applicable. The lu_decomposor should be used instead. 
		 */
		value_type det3 () const
		{
			const 	dynamic_vector<dynamic_vector<value_type> > & m = m_pElements ;
			value_type det = m[0][0] *(m[2][2]*m[1][1] - m[2][1]* m[1][2]) 
				- m[1][0] * (m[2][2]*m[0][1] - m[2][1]* m[0][2])
				+ m[2][0] * (m[1][2]*m[0][1] - m[1][1]* m[0][2]);
			return det ;
		}

		/**
		 * Inverse matrix of 3x3. 
		 */
		friend dynamic_matrix inv3 (const dynamic_matrix & m)
		{
			dynamic_matrix invm(3,3) ;
			float rdet = 1.0f / m.det3();
			invm[0][0] =   rdet * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) ;
			invm[0][1] = - rdet * (m[0][1] * m[2][2] - m[0][2] * m[2][1]) ;
			invm[0][2] =   rdet * (m[0][1] * m[1][2] - m[0][2] * m[1][1]) ;

			invm[1][0] = - rdet * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) ;
			invm[1][1] =   rdet * (m[0][0] * m[2][2] - m[0][2] * m[2][0]) ;
			invm[1][2] = - rdet * (m[0][0] * m[1][2] - m[0][2] * m[1][0]) ;

			invm[2][0] =   rdet * (m[1][0] * m[2][1] - m[1][1] * m[2][0]) ;
			invm[2][1] = - rdet * (m[0][0] * m[2][1] - m[0][1] * m[2][0]) ;
			invm[2][2] =   rdet * (m[0][0] * m[1][1] - m[0][1] * m[1][0]) ;
			return invm; 
		} 
	private:
		void allocate_memory ()
		{
			allocate_memory (m_pElements, m_row, m_col, value_type(0));
		}
		void allocate_memory (dynamic_vector<dynamic_vector<value_type> > & _m, unsigned _r, unsigned _c , const value_type & _v)
		{
			if (!(_m.size() == 0 ))
			{
				_m.clear();
			}
		
			dynamic_vector<value_type> temp(_c);
			for (unsigned i = 0; i < _r; ++i)
			{
				_m.append(temp);
			}
			for (unsigned i = 0; i < _r ; ++i)
			{
				for (unsigned j = 0; j < _c; ++j)
				{
					_m[i][j] = _v;
				}
			}
		}
		void deallocate_memory()
		{
			//this->deallocate_memory (this->m_pElements, m_row, m_col);
		}

	private:
		unsigned m_row, m_col ;

		dynamic_vector<dynamic_vector<value_type> > m_pElements ;


		/**
		 * For overloading comma. 
		 */
		unsigned m_idxcount ;  
	};
};




#endif
